Setup

First, load the libraries needed for the workflow

library(tidyverse)
library(CoRC)

Input

Put in all variables needed to customise the workflow

Put in the path to the model (SBML), as url or local path (as string) and load the model. Use loadModel() for .cps models.

modelPathSBML <- "Models/Leloup2004Circadian_model.xml"
model         <- loadSBML(modelPathSBML)

To use the results from the Sigma Point method, I will load the RDS with the results.

SP_result <- readRDS("docs/SP_Leloup")
SP_means  <- SP_result$sp_means

# set the Global Quantities to the estimated means
setGlobalQuantities(key = c("KAP", "KAC", "KIB"), 
                    initial_value = SP_means[c("Values[KAP].InitialValue",
                    "Values[KAC].InitialValue",
                    "Values[KAP].InitialValue")])

# read out Variance of Global Quantities
SP_cov    <- SP_result$sp_cov_matrix
SP_var    <- c()
for (i in 1:nrow(SP_cov)){
  SP_var[i] <- SP_cov[i,i]
}

Load the data into the variable data.

data <- read.csv("Data/Leloup2004Circadian_data.csv", sep = "")

Define experiment(s) for Parameter Estimation. For multiple experiments use list of defineExperiments().

pe_experiments <- defineExperiments(
  data = data,
  type = c("time", "dependent", "dependent", "dependent", "dependent", "dependent", "dependent"),
  mapping = c(NA, getSpeciesReferences(key = c("MP", "MC", "MB", "PTot", "CTot", "BTot"))$concentration),
  weight_method = "mean_square"
)

Set Parameter Estimation Parameters. For more parameters, copy the defineParameterEstimationParameter() function and append the list. All parameters in this model are mapped to global Quantities. For Parameter Estimation, the Global Quantities not assigned or with an ode (so only “fixed” global Quantities) are estimated.

globalQ <- getGlobalQuantities()[getGlobalQuantities()$type == "fixed",]$key

pe_parameters <- lapply(globalQ, make_pe_param_global, lower = 1, upper = 1)

Set Parameter Estimation Settings.

pe_randomize_start_values <- FALSE
pe_update_model           <- FALSE

pe_method <- setParameterEstimationSettings(
  method = "LevenbergMarquardt"
)

Set Parameter for Profile Likelihood. For more parameters, copy the defineParameterEstimationParameter() function and append the list, or use the make_pe_param() function defined in utils and run it on the parameters you want to add (see above). I want to use the variance as an upper and lower bound.

pl_interest   <- getGlobalQuantities(c("KAP", "KAC", "KIB"))$key
pl_parameters <- list(make_pe_param_global_absolute(pl_interest[3], 
                                                    upper = SP_var[3], 
                                                     lower = SP_var[3]),
                       make_pe_param_global_absolute(pl_interest[2],
                                                     upper = SP_var[2],
                                                     lower = SP_var[2]),
                       make_pe_param_global_absolute(pl_interest[3],
                                                     upper = SP_var[3],
                                                     lower = SP_var[3])
                      )

Choose maximum number of steps for Profile Likelihood.

max_steps     <- 50

Do you want to save the result as a RDS object? Also give the file, where you want to save it.

save <- TRUE
file <- "docs/PL_Leloup"

Workflow

This should generally not need to be changed, only if the workflow is to be changed more generally.

# Run Profile Likelihood 
PL_result <- runProfileLikelihood(
  model                     = model,
  pl_parameters             = pl_parameters,
  max_steps                 = max_steps,
  pe_parameters             = pe_parameters,
  pe_experiments            = pe_experiments,
  pe_method                 = pe_method,
  pe_update_model           = pe_update_model,
  pe_randomize_start_values = pe_randomize_start_values
)
#> [1] "Creating profile of parameter {Values[KIB].InitialValue}"
#> [1] "Creating profile of parameter {Values[KAC].InitialValue}"
#> [1] "Creating profile of parameter {Values[KIB].InitialValue}"

# Save result
if (save){
  saveRDS(PL_result, 
          file = file)
}

Output

This workflow produces the following output:

printLikelihoodProfile(PL_result)
#> [[1]]

#> 
#> [[2]]

#> 
#> [[3]]

Utils

Creates a function for easy Parameter Estimation Parameter Usage:

make_pe_param  <- function(name, upper, lower) {
  value_ref <- parameter_strict(name, reference = "Value")
  value     <- getValue(value_ref)
  defineParameterEstimationParameter( ref         = value_ref,
                                      start_value = value,
                                      lower_bound = abs(value)*lower,
                                      upper_bound = abs(value)*upper)
  
}

make_pe_param_global <- function(name, upper, lower){
  value     <- getGlobalQuantities(key = name)$initial_value
  value_ref <- getGlobalQuantityReferences(key = name)$initial_value
  defineParameterEstimationParameter(ref         = value_ref,
                                     start_value = value,
                                     lower_bound = value-abs(value)*lower,
                                     upper_bound = value+abs(value)*upper
                                     )
  
}

make_pe_param_global_absolute <- function(name, upper, lower){
  value     <- getGlobalQuantities(key = name)$initial_value
  value_ref <- getGlobalQuantityReferences(key = name)$initial_value
  defineParameterEstimationParameter(ref         = value_ref,
                                     start_value = value,
                                     lower_bound = value-lower,
                                     upper_bound = value+upper)
  }